Entomological parameters and population structure at a microgeographic scale of the main Colombian malaria vectors Anopheles albimanus and Anopheles nuneztovari

Population subdivision among several neotropical malaria vectors has been widely evaluated; however, few studies have analyzed population variation at a microgeographic scale, wherein local environmental variables may lead to population differentiation. The aim of the present study was to evaluate the genetic and geometric morphometric structure of Anopheles nuneztovari and Anopheles albimanus in endemic localities of northwestern Colombia. Genetic and phenetic structures were evaluated using microsatellites markers and wing geometric morphometrics, respectively. In addition, entomological indices of importance in transmission were calculated. Results showed that the main biting peaks of Anopheles nuneztovari were between 20:00 and 22:00, whereas Anopheles albimanus exhibited more variation in biting times among localities. Infection in An. nuneztovari by Plasmodium spp. (IR: 4.35%) and the annual entomological inoculation rate (30.31), indicated high vector exposure and local transmission risk. We did not detect Plasmodium-infected An. albimanus in this study. In general, low genetic and phenetic subdivision among the populations of both vectors was detected using a combination of phenotypic, genetic and environmental data. The results indicated high regional gene flow, although local environmental characteristics may be influencing the wing conformation differentiation and behavioral variation observed in An. albimanus. Furthermore, the population subdivision detected by microsatellite markers for both species by Bayesian genetic analysis provides a more accurate picture of the current genetic structure in comparison to previous studies. Finally, the biting behavior variation observed for both vectors among localities suggests the need for continuous malaria vector surveys covering the endemic region to implement the most effective integrated local control interventions.

Introduction Malaria remains one of the most important public health problems worldwide. The Plasmodium parasites that cause the disease are transmitted to humans by the female mosquitoes of approximately 40 species of Anopheles [1]. Colombia is third in the number of malaria cases in the Americas [2], and the number of cases/year in the past five years has fluctuated between 50,000-80,000 [3][4][5][6][7]. Currently, the most malaria endemic regions are the Pacific (PAC) in western Colombia (49.2%), and the Urabá-Bajo Cauca-Alto Sinú (UCS), in the northwest (18.7%) [3]. During the study period (2013)(2014), the number of malaria cases registered in UCS were 22,889 (2013) and 8,620 (2014), corresponding to 43.84% and 20.97% of the cases in the country, respectively [8,9]. The main malaria vectors in Colombia are Anopheles darlingi, Anopheles nuneztovari and Anopheles albimanus [10,11]. Vector control is one of the main strategies to decrease malaria incidence [1]; thus, knowledge of vector biology remains essential to reduce malaria transmission. It is known that anthropogenic environmental alterations and insecticide selection pressures affect vector population dynamics, because these factors may increase the abundance of anthropogenic Anopheles species and the appearance of insecticide resistant species, which affect malaria risk parameters [12][13][14]. Hence, genetic population structure and entomological parameters should be regularly monitored to evaluate changes in transmission risk to implement Integrated Vector Management programs [15]; such information will allow, for example, the identification of season(s) when vector control interventions should be intensified.
Population structure studies of the neotropical malaria vectors An. nuneztovari and An. albimanus have identified population subdivision across their distributions [16,17]. Regarding An. nuneztovari, several lineages were identified among South American populations [18], and a new species of the Nuneztovari Complex, lineage III, has been identified east of the Amazon region [19]. In Colombia, genetic differentiation and population subdivision was detected between An. nuneztovari populations from the west-northwest and the east-northeast, attributed to physical barriers, geographic distance and ecological variation on both sides of the Andes [20]. However, a study using COI and the white gene found high gene flow and the existence of a single taxon in An. nuneztovari from UCS, northwest Colombia [21]. Various studies support the status of An. albimanus as a single taxon in Central and South America [22,23], and genetic structure was identified between Central and South American populations [24]. In Colombia, An. albimanus is mainly distributed on the Atlantic and Pacific coasts [11,25], where the populations show little genetic differentiation [26,27]. The combination of genetic and morphometric data found low genetic structure and high gene flow, confirming An. albimanus as a metapopulation [26]. Both species have been detected infected with Plasmodium spp., An. nuneztovari in UCS-northwest and in PAC-west Colombia, and An. albimanus in PAC [25,28,29].
Previous studies in the Colombian malaria endemic UCS suggest that anthropogenic activities such as mining, livestock production, timber extraction, small scale agriculture and wood exploitation favor the presence of the main vectors An. albimanus and An. nuneztovari. These activities promote landscapes that impact key mosquito parameters, such as mosquito species composition, vector abundance, biting behavior and natural parasite infection, which determine malaria transmission risk at the local and regional levels [12,30,31]. The present study evaluates geometric morphometrics and genetic population variation in addition to entomological parameters for An. nuneztovari and An. albimanus at a microgeographic scale in the UCS region of northwest Colombia. These analyses may provide important information on the factors determining mosquito population variation that can influence vector transmission dynamics.

Study area and sample collection
Anopheles specimens were collected from two localities in each of five municipalities, six sites per locality, in the malaria endemic region Urabá-Bajo Cauca and Alto Sinú (UCS) in northwest Colombia (Table 1, Fig 1). Selection of the localities and collection times was based on reports of malaria transmission and safety considerations. Collections were performed from September 2013 to October 2014, during three nights per locality from 18:00-24:00 h, in open areas of the house and in the peri-domestic area using protected humans as an attractant (Informed consent agreement and collection protocol reviewed and approved by a University of Antioquia Institutional Review Board -Bioethics Committee, Facultad Nacional de Salud Pública-Universidad de Antioquia, Acta 063). A written informed consent was obtained from all individuals participating in the collections. The specimens were identified using a morphological key [10] and species assignment was confirmed by PCR-RFLP-ITS2 [32][33][34] and COI barcode gene sequencing [35,36].

Analysis of genetic data
Given the geographical proximity among localities from the same municipality, specimens were grouped and analyzed by municipality for the genetic analysis. Inter and intra-population genetic diversity was estimated for 139 An. nuneztovari specimens and 141 An. albimanus (Table 1). For An. nuneztovari, 11 microsatellite (MS) loci were analyzed (Anu3, Anu7, Anu11, Anu12, Anu10, Anu1, Anu2, Anu9, Anu4, Anu6, Anu8) [19], (S1 Table) [40]. Genetic population structure at the population and individual levels was analyzed using a grouping method based on both transient Hardy-Weinberg disequilibrium (HWD) and linkage disequilibrium (LD) caused by admixture between populations, in Structure v.2.3.4 [41]. Based on allele diversity data, individuals with unique alleles were grouped together into assumed populations (K) which is pre-determined [42]. The K value with the maximum posterior probability Pr (X|K) was retained and assumed to be the most probable number of clusters in that population. Twenty independent runs were performed for each value of K (K = 1 to 10) with a burnin period of 500,000 iterations and 2 million steps for Monte Carlo Markov Chain (MCMC) replications run using Structure Harvester v.0.56.3 [43]. A total of 25 files with optimal K were averaged using the Large K Greedy method implemented in Clumpp v. 1.1.2 [44]. A Cavalli-Sforza chord distance and Neighbor Joining (NJ) consensus tree representing genetic differentiation among populations was created using Populations v. 1.2.31 [45] with 1000 bootstrap replicates [46]. Population differences defined by Structure were visualized using a Factorial Correspondence Analysis (FCA) in Genetix v. 4.05.2 [47]. The number of alleles (Na), expected heterozygosity (He), observed heterozygosity (Ho), allele richness (Rs), number of migrants per generation (Nm), inbreeding index (F IS ), fixation index (F ST ) and balance of Hardy-Weinberg (HWE) were calculated for each MS loci using Arlequin software v. 3.5 [48]. The number of private alleles for all loci was obtained using Convert v. 1.31 [39]. Molecular analysis of variance (AMOVA) was used to examine within and among populations variation in Arlequin; the allelic frequencies were used to estimate Euclidean distances at different levels of genetic structure, with 10,000 non-parametric permutations [48].

Morphometric data analysis
Wing size and conformation variation were evaluated by geometric morphometric analysis. The right wings of 147 An. nuneztovari specimens (27-30 per population) and 146 An. albimanus (26-30 per population) ( Table 1), were mounted on microscope slides with coverslips and photographed using a digital camera (Moticam1 2500), attached to an Olympus 1 Sz61 stereomicroscope. Right wing images were used for morphometric analysis and a set of 13 reference measurements were considered in the analysis [26]. To estimate variation of wing size and conformation, raw data coordinates were overlaid using a Procrustes full fitting procedure that eliminates variation due to scale, position and orientation of datum settings [49]. Wing size variation was examined using the centroid size (CS), defined as the square root of the sum of the square distances of a set of reference points from the centroid [50]. CS variation was tested using a Mann-Whitney test after Bonferroni correction [51].
Wing conformation among populations was compared by means of a disparity analysis, using 10,000 random permutation tests. To calculate the differences in wing conformation among populations, allometry-free variables were used as input for a Canonical Variates Analysis (CVA). To obtain Procrustes and Mahalanobis distances, the CVA analysis was used to determine whether the five geographic populations of each species could be statistically distinguished based on the relative deformation matrix [26,52]. Means of the wing conformation for specimens of the five geographic populations were plotted along the two axes of canonical variation based on the Procrustes distance matrix [26,52]. Reference point digitization, morphometric analysis and graphical results were performed using various modules of the Clic package [53]. The Past program was used for wing size comparisons [54].

Relationship among phenotypic, genotypic and environmental data
Population structure based on allele frequency variation of georeferenced individual genotypes and inference of the number of populations and spatial location of genetic discontinuities between those populations was evaluated by a Bayesian clustering algorithm in Geneland 4.0.3 [55], implemented in R software v. 3.3.2 (R Development Core Team 2008). To estimate the optimal number of subpopulations, based on the spatial location of the sampling sites, a set of data was analyzed using combinations of the phenotypic (P), genetic (G) and spatial (S) matrices, as follows: 1) phenotypic data under the spatial model; 2) genetic data under the spatial model; and 3) phenotypic and genetic data under the spatial model [24]. Potential spatial distribution models (Spatial distribution model-SDM) were generated in Maxent v. 3.3.3 [56], using as input the layers of the Normalized Difference Vegetation Index (NDVI) and An. nuneztovari and An. albimanus occurrence records throughout the endemic region, to generate an environmental resistance matrix in Circuitscape v. 3.5.8 [57]. Finally, the relationship between environmental heterogeneity with geographic and with genetic distances, and phenotypic variation of the population, was evaluated using the partial Mantel test [58]. The matrices of genetic structure (F ST ), phenotypic differentiation (Mahalanobis distances), environmental or cost distances and geographical distances were compared using R v. 3.3.2.

Entomological parameters and entomological inoculation rate for Plasmodium
Mosquito parasite infection was tested in the pools of five mosquito heads and thoraxes by an enzyme-linked immunosorbent assay (ELISA) using monoclonal antibodies against Plasmodium falciparum, Plasmodium vivax VK247 and VK210. The positive pools were confirmed by a second ELISA and a nested PCR performed with DNA extracted from the mosquito abdomen and Plasmodium specific primers [29,59]. The entomological parameters were estimated by locality and included, biting behavior that was calculated as the number of mosquitoes of each species collected by hour and locality. Human Biting Rate (HBR), estimated as the average number of Anopheles females collected per hour and per average number of collectors, expressed as the number of bites per person, per night (b.p.n) [60]. Infection rate (IR) was calculated as the percentage of Plasmodium positive mosquitoes out of the total number of mosquitoes analyzed, by species and locality. The annual entomological inoculation rate (EIR) per site that corresponds to the number of infective bites that a person may receive in one year [29].

Genetic diversity and population genetic structure
Genetic diversity analyses for An. nuneztovari showed that the ANU11 locus was not polymorphic and was thus excluded from population analyses. For An. albimanus all loci were polymorphic. In An. nuneztovari, the expected mean loci heterozygosity (He) ranged from 0.613 in El Bagre to 0.674 in Mutatá and for An. albimanus from 0.802 in Turbo to 0.840 in Arboletes (Table 2); also, the number of alleles, observed heterozygosity and inbreeding coefficient differs between both species but not as much between localities for each species ( Table 2). The F ST and Nm values indicated low genetic differentiation and high gene flow among An. nuneztovari populations in the UCS region. The average Nm was high, with the highest value between Cáceres and Tierralta (64.5) (S3 Table). Bayesian clustering analysis in Structure revealed three subpopulations of genetically related individuals (K = 3). There was no particular association among these subpopulations and geographic location (Fig 2A). Similarly, F ST and Nm values for An. albimanus in UCS showed low genetic differentiation and high gene flow among populations. The average Nm was high, with the highest value detected between Moñitos and San Antero (526.4) (S4 Table). The genetic structure analysis revealed two subpopulations of genetically related individuals (K = 2); however, there was no association among these subpopulations and geographic location (Fig 2B).

Population phenotypic variation
Geometric morphometry analysis of 60 An. nuneztovari wings showed good repeatability in the x, y coordinates (R = 0.953), centroid size (R = 0.899) and relative deformations   Table). Paired comparisons using the Mann-Whitney test showed significant differences for various centroid sizes (Fig 3B, S6 Table), mainly for the Arboletes population compared to the four others.
The metric disparity analysis for An. nuneztovari wing conformation showed significant differences between the following comparisons: Cáceres and Mutatá (p = 0.0145), Cáceres vs El Bagre (p = 0.0007), Tierralta and Mutatá (p = 0.0433), Tierralta and El Bagre (p = 0.0140) and Turbo and El Bagre (p = 0.0152). Discrimination among populations in the morphospace showed that the first two CV axes represented 76% of the total variation of the data (Fig 4AI). Genetic differences among populations analyzed by FCA showed that axes 1 and 2 explained 39.07 and 24.49%, respectively ( Fig 4AII); this analysis showed a slight separation of the El Bagre population. For An. albimanus the wing shape metric disparity analysis revealed no significant differences between the paired comparisons. The discriminant analysis among populations in the morphospace showed that the first two CV axes represented 75% of the total variation of the data (Fig 4BI). The genetic differences among An. albimanus populations analyzed by FCA showed that axes 1 and 2 explained 2.70 and 2.68% respectively (Fig 4BII). The phenotypic discriminate analysis showed a slight separation of the Arboletes population, observed to a lesser degree in the genetic analysis.

Combination of phenotypic, genetic and environmental data
Bayesian analysis of genetic data in GENELAND, which considers the spatial distribution model, showed the presence of five An. nuneztovari genetic subpopulations (Fig 5AI). However, the same analysis with only the phenotypic data showed that there are two subpopulations ( Fig 5AII); the phenotypic-genetic data corroborated the existence of five subpopulations but with a different geographic distribution (Fig 5AIII). A significant correlation was found between genetic vs. geographic distances by the Mantel test (S7 Table). Three An. albimanus genetic subpopulations were identified by Bayesian analysis (Fig 5BI). However, the same analysis with only phenotypic data showed only one subpopulation (Fig 5BII), while the phenotypic and genetic data revealed two subpopulations (Fig 5BIII). A Mantel test showed no significant correlations among any of the distance matrices (S8 Table).

Biting behavior
Intra-and interspecific differences were observed in An. nuneztovari (Fig 6) and An. albimanus biting behavior (Fig 7). For An. nuneztovari, the main biting peaks were between 20:00h and 22:00 h, however, differences were observed in the Montelibano localities with a main peak at 18:00 h in Puerto Nuevo and at 21:00 h in Puerto Anchica. A late peak at 23:00 h was observed in Asturias, Cáceres municipality (Fig 6). The localities Camerún in Turbo and Rio Cedro in Moñitos were not included in the biting behavior analysis because of low mosquito density. Variation in biting times was observed for An. albimanus among the localities, however, two main biting peaks were evident: one in the early evening (18:00-20:00) h in Bahia Cispatá and Puerto Nuevo and another from 20:00-23:00 h in Camerún-Turbo, Broqueles-Moñitos, La Arenosa and the main settlement of Arboletes (Fig 7).

Entomological parameters
The human biting rates (HBR) in most localities were relatively homogeneous with rates up to 6 bites per person per night (b.p.n) for An. nuneztovari and 9 b.p.n. for An. albimanus ( Table 3). The highest HBRs for both species were registered in Asturias-Caceres for An. nuneztovari (5.91 b.p.n) and in La Arenosa-Arboletes for An. albimanus (53.6 b.p.n; Table 3). The lowest rates were detected in Camerún-Turbo (0.25 b.p.n.) for An. nuneztovari and Santa Ana-Tierralta (0.33 b.p.n.) for An. albimanus. The infection rate was of 4.35%, given the detection of one An. nuneztovari specimen infected by Plasmodium spp., in Puerto Anchica-Montelibano, for an annual EIR of 30.31 infective bites per year (Table 3).

Discussion
This study evaluated genetic population structure, geometric morphometrics and relevant entomological parameters for An. nuneztovari and An. albimanus at a microgeographic scale in the UCS malaria endemic region. A previous study of An. nuneztovari in Colombia showed that this species is unique, with a population subdivision between the west-northwest and the east-northeast, which suggested high genetic differentiation and reduction of gene flow among populations caused by distance, ecological differences and the presence of geographic barriers (i.e., the Andes) [20]. In the present study, a Bayesian genetic analysis using microsatellite markers, which exhibit a high mutation rate, allowed for the identification of population variation at a microgeographic scale. Of note, three co-occurrent An. nuneztovari subpopulations were detected that presented low genetic differentiation and high gene flow. In contrast, previous studies of An. nuneztovari in UCS using mitochondrial COI and nuclear white gene markers suggested genetic homogeneity in northwestern Colombia [21]. The present study used microsatellite markers which are characterized as codominant, highly polymorphic, having Mendelian inheritance, and are easily typed making them suitable genetic markers for analyzing population structure, genomic variation and evolutionary processes [61]. The present results reflect the utility of the microsatellite markers in detecting population variation in An. nuneztovari, showing higher sensitivity compared with mitochondrial and other nuclear markers. Similarly, a previous study on An. nuneztovari s.l. from the Brazilian Amazon region based on the same microsatellite markers, reported the presence of three genetic lineages [19]. Given our results in comparison with earlier data, we suggest that microsatellite markers are providing a more accurate picture of the current population structure of An. nuneztovari in UCS.
Anopheles nuneztovari geometric morphometric analysis revealed significant differences in size and wing conformation for several of the paired comparisons between populations. This species has shown high wing morphological variability previously in specimens from northwest Colombia [62,63], Venezuela [64] and other regions of Latin America [65]. We found no relationship between the phenotypic and genetic distance matrix which indicates that the wing size differences detected among An. nuneztovari specimens are more likely associated with  phenotypic plasticity rather than wing conformation variation as previously suggested for other anophelines [66]. It is also possible that wing size differences among populations are related to altitudinal differences among populations (0-76 m.a.s.l.). A previous study on Anopheles cruzi form São Paulo state in Brazil indicated vertical population structuring of wing geometry despite the similarity among landscape microenvironments; wing shapes were distinct between lowland and hilltop populations and the wings of hilltop specimens were larger [67]. For An. albimanus from UCS, Bayesian analysis detected two subpopulations with low genetic differentiation, high gene flow and no evidence of geographic differences. The absence of any obvious physical and biogeographic barriers in the UCS region may explain the low genetic differentiation recorded. Previous studies at a larger scale showed moderate to low genetic differentiation between Colombian Caribbean and Pacific populations of An. albimanus and suggested that this metapopulation was influenced by demographic events during the late Pleistocene [26,27]. Furthermore, a panmictic population was detected between Central and South America [16,24]. The Gómez et al [26] and Gutiérrez et al., [27] studies of Colombian An. albimanus were conducted with four microsatellite markers; in this study 11 microsatellites were evaluated and provided a finer resolution as shown by the Bayesian analysis that was able to identify three genetic subpopulations. We hypothesize that the low genetic differentiation and high gene flow in the present study is influenced by the dispersion capacity of An. albimanus adults, with an average flight range of 32 km [68]; in addition, the absence of a notable physical barrier such as high mountains or wide rivers in the endemic area likely favored gene flow. We hypothesize that adult mosquito movement by human or wind-assisted transport favors displacement of An. albimanus to nearby populations increasing gene flow and decreasing population structure.
Of note, significant variation in centroid size was observed in all paired comparisons for An. albimanus from Arboletes in the geometric morphometric analysis. Also, the multivariate alar conformation and genetic analysis showed a slight separation of the Arboletes population. A previous study of An. albimanus populations from the Colombia Caribbean and Pacific Coasts suggested size variation occurred as a consequence of environmental change experienced during the immature stages [26]. We hypothesize that in our study, Arboletes larval environments are distinctive in some essential parameters for anopheline development such as water temperature, pH and/or oxygen availability, as detected by Soleimani-Ahmadi et al. [69] in their study of anopheline larval habitat differences in southern Iran. In addition, factors such as food availability, larval population density and predators may influence intrapopulation adult size variation of Diptera [70].
In addition to genetic and phenetic evaluations, entomological indices were calculated because these estimates are fundamental tools for the design of malaria vector control strategies [71]; they provide information about changes in vector population dynamics that affect malaria transmission risk, which is important for determining when and where control measures are necessary. Variation in biting behavior has been reported for An. nuneztovari and An. albimanus across their Latin American distributions [28,29,72]. In this work, biting peaks for An. nuneztovari varied according to the locality; however, its main biting peaks in most localities were between 20:00-22:00 h. In other localities of this region, a similar biting peak for An. nuneztovari was observed [28, 30,73]. Altogether, our results suggest the need to conduct regular surveys to monitor this vector behavior to appropriately direct control measures. For An. albimanus, biting peaks also varied, but the highest was found in Arboletes between 22:00-23:00 h. These values differ from those previously reported, in 2009, for this endemic area (An. albimanus 18:00-22:00 h) [73]; perhaps some environment parameters have been altered as a result of increasing anthropogenic activity in the region, influencing the biting time of this species. Variation in An. nuneztovari and An. albimanus biting behavior was not associated with genetic population differences, and several studies have demonstrated the effects of plasticity on Anopheles species biting behavior [74,75]. Similar to the results of the present study, a recent investigation of An. darlingi from the Amazon region found an association of genetic diversity with biting behavior not related to population structure [76]. In addition, the variation in biting peaks for both vector species may be related with the degree of human exposure [28]. In this context, humans are under higher risk of being bitten by Plasmodium infected Anopheles when they carry out leisure activities in the open spaces of their houses, at the hours when the vectors show their higher biting activity.
In conclusion, results of entomological parameter estimations indicated that An. nuneztovari and An. albimanus exhibited variable biting behavior among localities and suggest that the UCS human populations are exposed continuously to malaria vector bites mainly during the early part of the night; therefore, control measures should focus to reduce vector-human contact at this time. In addition, the Plasmodium infection results for An. nuneztovari indicate that even the low human biting rate (< 2 b.p.n.) is sufficient to maintain malaria transmission, and although, An. albimanus is an important malaria vector, it was not detected infected; this can be due to the low infection rates usually detected in endemic regions of Colombia [28]. Furthermore, the low population structure detected for both vectors indicates high gene flow in the region, although local environmental characteristics may influence the wing conformation differentiation and behavioral variation. Our results indicate the importance of evaluating population structure and behavior at the local level to design the most cost-effective, targeted control strategies.
Supporting information S1